3D projection method

ABSTRACT

A digitized tomosynthesis method for projecting a 3D volumetric image of an object onto a virtual projection plane that comprises an alpha blend plane subdivided into a rectangular array of pixel cells. A ray of energy defines a trace representing a virtual sight path ray passing from a view point through one pixel cell of the alpha blend plane to define a view level image, wherein a 3D volumetric image of the object is projected by computing the coordinate of an object voxel referenced to the image plane, computing the image plane intercept of the ray of energy from the source through the object voxel, and forming a plurality of successively selected view level image planes at the alpha blend plane.

CROSS REFERENCE TO RELATED APPLICATION

This application claims the benefit of Provisional Patent ApplicationNo. 60/422,532, filed Oct. 31, 2002, and is a continuation-in part ofU.S. patent application Ser. No. 09/709,586 filed on Nov. 13, 2000 nowU.S. Pat. No. 6,671,349, as well as being a continuation-in-part of U.S.patent application Ser. No. 10/687,180 filed on Oct. 15, 2003, which inturn is a continuation-in part of U.S. patent application Ser. No.09/709,586 filed on Nov. 13, 2000.

FIELD OF THE INVENTION

The fields of art to which the invention pertains include the fields ofdynamic tomography and computed tomography.

BACKGROUND OF THE INVENTION

A unique x-ray imaging and inspection system that combines 3D volumetricimaging and conventional 2D radiography for a complete x-ray inspectionsolution has been provided under the mark Digitome. For purposes of thisapplication, the system will be referred to a “digitized tomography”.Digitized tomography technology has been used for film based3-dimensional x-ray imaging, and has been enhanced to incorporatedigital flat panel x-ray detectors, resulting in examinations being madein minutes. Its features, provide unique capabilities to view anyhorizontal or vertical plane, scan through the volume in 0.005″increments and measure internal features. See Griffith U.S. Pat. No.5,051,904 (“Computerized Dynamic Tomography System”), U.S. Pat. No.5,070,454 (“Reference Marker Orientation System For A RadiographicFilm-Based Computerized Tomography System”), and U.S. Pat. No. 5,319,550(“High Resolution Digital Image Registration”), the disclosures of whichare incorporated herein by reference.

The digitized tomography software contains what is known as thedigitized tomography kernel. It is a set of software modules that areused to compute digitized tomography views. One defines the geometry bywhich the images are formed, the data arrays that contain the imagedata, and the coordinates of an object space. The resultant voxel valueis returned to the calling software. It is the calling software'sresponsibility to resend the coordinates and store the voxel values sothey produce the desired digitized tomography view foil.

From its inception, the digitized tomography kernel was based upon adigital simulation of what is known as the film mode. Referring to FIG.1, its geometric derivation assumes the source 10 is at a canted anglewith respect to the image plane 12. The object to examine is positionedcentrally over the intersect 14 of the x-ray source optical axis 16 andthe image plane. The object is then is stepwise rotated about an axisperpendicular to the image plane 12 and that passes through the aboveintersect. An x-ray image is acquired at each rotation position of theobject.

In the film simulation mode, each acquired image is rotated by the sameangle as was the object when it was exposed. The resulting images arestacked with their axes of rotation coincident. Finally, each image istranslated radially inward from its axis of rotation by a distancerelated to the object level one desires to observe. The translationdistance is a function of the perpendicular distance between the imageplane and the view plane. Light passing through a stack of filmpositioned thusly will reveal the features at the desired view plane.

A similar computed combination of the data from a series of digitallyacquired images will result in the same view obtained in the previouslyreferred to film simulation mode and requires the combination ofregistered pixels from each of the digitized images. The method isdescribed in Appendix A. In addition, other view foils are possible tocompute simply by changing the pattern of the voxel coordinates passedto the digitized tomography kernel as well as various imageenhancements. The film mode simulation digitized tomography kernel mathis described in Appendix B wherein the kernel is limited to having theaxis of rotation of the object perpendicular to the image plane.

Referring to FIG. 2, an alternate geometry can be useful. For example,the x-ray optical axis 16 could be perpendicular to the image plane andthe axis of rotation of the object could be at another angle withrespect to the image plane 12. This can allow the instrumentation to fitinto a much smaller column because the source is immediately above theimage plate rather than offset to the side. It can also permit largerand taller objects to be examined using a given image plate than doesthe current geometry. In a film mode simulation geometry, the objectplane is not parallel to the film plane and the geometric magnificationvaries over the film plane. Methods can be provided by which a markerdisk and misalignment post are placed to calibrate the alternategeometry.

BRIEF SUMMARY OF THE INVENTION

The present invention is directed to a 3D projection method for adigitized tomography image stack that includes the geometry of FIG. 1,the alternate geometry of FIG. 2, and potentially many other geometries,and uses what can be referred to as the digitized tomography kernel. Thepresent invention uses the digitized tomography image stack to generatea 3D projection, using methods described more specifically in theDetailed Description below and in appendices A, B, C and D.

More particularly, a digitized tomosynthesis method is provided forprojecting a 3D volumetric image of an object onto a virtual projectionplane. The image is defined by object voxels, in which a ray of energyfrom a source, preferably electromagnetic radiation such as x-rayradiation, travels through the object to impinge on an energy sensingimage plate, such as a flat panel digital detector, having an activearea defining an image plane and in which an image is acquired by theenergy sensor at successive relative rotational positions of the objectand source. The image plate is parallel with a defined XZ plane and thesource is positioned above the XZ plane with its optical axis definingthe center of a radiation cone and intersecting the image plate atapproximately its central pixel. A 3D volumetric image of the object isprojected by computing the coordinate of an object voxel referenced tothe image plane, computing the image plane intercept of the ray ofenergy from the source through the object voxel, and projecting areconstructed 3D volumetric image of the object onto the virtualprojection plane.

In a specific embodiment, the coordinate of the object voxel referencedto the image plane is computed by (a) defining an object bounding boxabove the XZ plane that intercepts the radiation cone so that a shadowof the object falls on the active area of the image plate, the positionand orientation of the bounding box having a reference point coordinatedefined by X, Y, Z, pitch, yaw, and roll, (b) specifying the coordinateof the voxel referenced to the object bounding box, (c) rotating thevoxel coordinate by the pitch, yaw, and roll of the bounding box, and(d) translating the voxel coordinate to the reference point coordinate.The bounding box has a size and height whereby it includes the volume ofinterest of the object. The bounding box and the image plate can each berectangular, the sides of the bounding box being parallel with the sidesof the image plate.

The virtual projection plane can comprise an alpha blend planesubdivided into an array, e.g., a rectangular array, of pixel cellswherein the ray of energy defines a trace thereof representing a virtualsight path ray passing from a view point through one pixel cell of thealpha blend plane to define a view level image. There is a ray trace foreach pixel cell of the alpha blend plane for each view level image. Aprojected planar view of a selected view level image plane is formed atthe alpha blend plane by transferring the image pixel cell value at theintersect of the ray trace and the selected view level plane to theimage pixel cell at the intersect of the ray trace and the alpha blendplane. A plurality of view level images are formed successively at oneimage plane, the view point and alpha blend plane being translated alevel closer to said one image plane for each successive view level.

The 3D volumetric image can be stored in a computer memory as an arrayof pixels represented by a value in which the brightness of each pixelis proportional to the value. The steps of the method of this inventioncan be performed mathematically and programmed whereby the reconstructed3D volumetric image is projected by execution of the program on acomputer.

Features and advantages of the invention will be described hereinafterwhich form the subject of the invention. It should be appreciated bythose skilled in the art that the conception and specific embodimentdisclosed may be readily utilized as a basis for modifying or designingother structures for carrying out the same purposes of the presentinvention. It should also be realized by those skilled in the art thatsuch equivalent constructions do not depart from the spirit and scope ofthe invention. The novel features which are believed to becharacteristic of the invention, both as to its organization and methodof operation, together with further objects and advantages will bebetter understood from the following description when considered inconnection with the accompanying figures. It is to be expresslyunderstood, however, that each of the figures is provided for thepurpose of illustration and description only and is not intended as adefinition of the limits of the present invention.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic depiction of film mode operation of the digitizedtomography kernel wherein the source is at a canted angle with respectto the image plane;

FIG. 2 is a schematic depiction of an alternate geometry for operationof the digitized tomography kernel wherein the x-ray optical axis isperpendicular to the image plane and the axis of rotation of the objectis at another angle with respect to the image plane;

FIG. 3 is a schematic depiction of a method of the present invention forprojecting a 3D image;

FIG. 4 is a schematic representation of a film mode simulation geometrywherein the object plane is parallel to the film plane and the geometricmagnification projects uniform shadows of the features of the objectplane;

FIG. 5 is a schematic representation of an alternative geometry whereinthe object plane is not parallel to the film plane and the geometricmagnification varies over the film plane;

FIG. 6 depicts a marker in the existing system in which the object is around marker disk with a central pin made of a material suitably opaqueto x-rays;

FIG. 7 depicts use of the disk and pin of FIG. 6 in a first step incalibrating the geometry of the system;

FIG. 8 depicts a plastic post with a similar pin at one end and hole toreceive the pin of the round disk, inserted into the central hole of theobject turntable in a second step in calibrating the geometry of thesystem;

FIG. 9 depicts one method by which the marker disk and misalignment postare placed to calibrate the alternate geometry;

FIG. 10 depicts another method by which the marker disk and misalignmentpost are placed to calibrate the alternate geometry;

FIG. 11 shows the representation of 3D space in a “left-handed”coordinate system, used in 3D mathematics;

FIG. 12 illustrates a plane showing normal as the direction that thesurface of the plane is facing, used in 3D mathematics;

FIG. 13 illustrates the calculation of a plane from 3 given points, usedin 3D mathematics;

FIG. 14 illustrates 3 points stored in a clockwise direction in the x/yplane, used in 3D mathematics; and

FIG. 15 illustrates mirroring wherein all points are flipped about anarbitrary plane in 3D space, used in matrix mathematics.

DETAILED DESCRIPTION OF THE INVENTION

Referring to FIG. 3, the digitized tomography view level image planestack 11 represents a sequence of digitized tomography planar viewsgenerated by methods disclosed in the appendices A, B, C and D. Theimages 19 on the stack 11 are generated separated by a fixeddistance—usually by one pixel. The axis of generation can be anyarbitrary axis. The most advantages axis is coincident with the examinedobject's axis of rotation otherwise known as the “Z” axis. For such asituation the planes can be generated as “XY” views. A view level imageplane stack 11 can contain the least number of views to represent theobject being projected. The “Z” and “XY” designations are due tohistorical decisions. The adaptation of modern graphical convention andray tracing methods necessitated a change of axis labels where Z becomesY, X becomes Z, and Y becomes X.

More specifically, the following defines a generalized method for thecomputation of the coordinate of the shadow of a specified digitizedtomography object voxel:

The Coordinate System

The coordinate system follows the left hand rule using an aircraft yaw,pitch, roll metaphor to describe rotations. The eye is looking in thepositive direction of the Z axis. The X axis is horizontal and positiveto the right. The Y axis is vertical and positive up. All positiverotations will be counter clockwise when viewed along the negative axisdirection. Pitch will be a rotation about the X axis. Yaw will be arotation about the Y axis. Roll will be a rotation about the Z axis.Normal rotation transformations apply rotations in pitch, yaw, roil (X,Y, Z) order.

The Coordinate System Applied to Digitized Tomography

Three things are accurately located in space: the image plate, theobject, and the source. The coordinates are to be specified in pixelunits. One pixel unit is the distance between pixels on rows or columnsof the image plate. The first pixel of the image plate is coincidentwith the coordinate axis origin and the plate is parallel with the XZplane. The source is positioned above the XZ plane and its optical axis(center ray of radiation cone) intersects the image plate atapproximately the image plate's central pixel. The radiation cone ispresumed to be of approximately uniform intensity and emitted as if froma point source.

The object bounding box is positioned above the XZ plane and interceptsthe radiation cone so that the shadow of the object falls on the activearea of the image plate. The bounding box base is a square. The size ofthe base and the height of the bounding box is to be such that itincludes the volume of interest of the object. Positions within theobject are referenced to the center of the base of the bounding box. Thestandard position of the bounding box is with its base plane centered onorigin of the image plate plane and parallel with the image plate plane.The sides of the bounding box are parallel with the sides of the imageplate. The final position and orientation of the bounding box isexpressed as the six degrees of freedom coordinates of the referencepoint: X, Y, Z, pitch, yaw, roll

Computing the Image Coordinate of a Digitized Tomography Object VoxelShadow

-   1. Compute the coordinate of the object voxel referenced to the    image plane    -   a. Specify the coordinate of the voxel referenced to the object        bounding box    -   b. Rotate the voxel coordinate by the specified pitch, yaw, and        roll of the bounding box    -   c. Translate the voxel coordinate to the reference point        coordinate.-   2. Compute the image plane intercept of the ray from the source    through the object voxel    3D Projection Method

The view point 13 represents the virtual position of the eye of theobserver. The alpha blend plane 15 represents the virtual projectionplane for the reconstructed projected 3D image. The alpha blend plane 15is subdivided into a rectangular array of cells (i.e., pixels). The raytrace 17 represents a virtual sight path ray passing from the view point13 through one pixel of the alpha blend plane 15 to one of the viewlevel images. There will be a ray trace 17 for each cell of the alphablend plane 15 for each view level image plane 19.

A projected planer view of any given view level image plane 19 can beformed at the alpha blend plane 15 by transferring the image cell valueat the intersect of the ray trace 17 and the given view level plane 19to the image cell at the intersect of the ray trace 17 and the alphablend plane 15. A projected 3D image on the alpha blend plane 15 can beformed by alpha blending successive projected planer views, formed asabove, starting from the bottom of the stack and starting with a blackalpha blend image plane.

The above process can be accomplished mathematically and programmed forexecution on a computer. However, a much more computationally efficientmethod is to form the images of the view levels successively at oneimage plane. Then translate the view point and the alpha blend plane onelevel closer to that plane for each successive view level.

In an implementation of the 3D projection method, an image is stored incomputer memory as a rectangular array of pixels. A pixel is stored asan integer with a sufficient number of bits to represent the dynamicrange of the image. The size of the pixel is determined by the imageacquisition device. The pixels are stored in increasingly higher memorylocations from left to right and bottom to top of the image. The bottomto top arrangement is a consequence of the way current operating systemshandle image arrays and is not a property of digital images as such. Thebrightness (whiteness) of a pixel is proportional to the value stored inthe pixel. White is represented as the maximum pixel value and black aszero. In the case of the current implementation, a pixel is stored as asixteen bit word with a range of 0 to 4095 (12 bits) and has a size of0.005 inches.

In the process of the present invention, the following steps areimplemented:

-   1. Form a view level plane by method described in Appendix E. The    origin of the plane corresponds with the first pixel of the image.    Increasing X axis corresponds with image bottom to top. The Y axis    is normal (perpendicular) to the plane. The increasing Z axis    corresponds with image left to right. A pixel position corresponds    to a integer Z, X coordinate-   2. Compute reference coordinate    -   Xr=image height/2    -   Yr=0.0    -   Zr=image width/2.-   3. Compute the coordinate of the view point from operator provided    data.    -   Xvp=Xr    -   Yvp=view distance/pixel size    -   Zvp=Yr    -   Rotate Xvp,Yvp,Zvp about the Z axes through Xr,Yr,Zr by (90        degrees—view angle) by method described in Appendix F.-   4. Set the view level 19 to start level-   5. Set all cells of an alpha blend plane 15 image data matrix to    black.-   6 For each view level until a view level greater than end level:    -   Yvp=(view distance−(start level−view level))/pixel size    -   For each pixel of the alpha blend plane        -   Xpp=image line index        -   Ypp=(projection plane distance−(start level−view            level))/pixel size        -   Zpp=image column index-   Rotate Xpp, Ypp, Zpp about the Z axes through Xr, Yr, Zr by (90    degrees—view angle) by the method described in Appendix F.    -   -   Compute intersect of ray through Xvp,Yvp,Zvp and Xpp,Ypp,Zpp            and the view level plane by method described in Appendix E.        -   Extract the image cell value from the view level image at            the intercept coordinates by method described in Appendix B.        -   Save extracted image cell value in temporary projected image            data matrix        -   a the cell designated by the alpha blend plane image line            index and image column index.

    -   Alpha blend the temporary projected image data matrix as source        with the alpha blend plane image data matrix by the following        method:

    -   Set view level to view level+view level increment.        Alpha Blendinq

Alpha blending is a standard graphics method for combining digitalimages or computer graphical elements such that there are varyingdegrees of transparency. The goal is to make background features visiblethrough selected foreground features while still being able to see theforeground feature. The calculation:destination pixel=(alpha*source pixel)+((1−alpha)*destination pixel)is applied to each pixel position of a paired source and destinationimage or graphic. An alpha value of 1.0 produces the effect of thesource pixel being opaque. An alpha value of 0.0 produces the effect ofthe source pixel being totally transparent. Alpha values between 1.0 and0.0 will produce the effect of the source pixel being translucentinversely proportional to alpha.

In “true” color systems that implement alpha blending, it is typical tohave 32 bit pixels that contain a red, blue, green, and alpha value foreach pixel. The four values are each contained in a byte and can assumevalues from 0 to 255. For the three color values, 0 represents black,255 represents the given color being of full brightness, and values inbetween represent less than full brightness proportional to its value.For the alpha value, 0 signifies the pixel is transparent and 255signifies the pixel is opaque. Alpha values between 0 and 255 representvarious levels of translucency inversely proportional to the alphavalue.

Digitized tomography uses monochrome pixels. In the current system, apixel is contained in sixteen bits (two bytes) and holds a twelve bitvalue from 0 to 4095. Where the 0 level represents a fully exposed pixeltypically presented as black. The 4095 level represents a totallyunexposed pixel typically presented as white. Since the image is theresult of an X-ray exposure, the higher the pixel value the more dense(X-ray opaque) the object being exposed.

To blend X-ray images so the more dense features show through the lessdense features, black pixels must be transparent and white pixels mustbe opaque. Grey must be treated increasingly translucent the blacker itbecomes. To accomplish this effect, alpha is computed as a function ofthe source pixel value:alpha=(source pixel/(maximum pixel value+1)).In current systems the maximum pixel value is 4095. A configurableexponent is added to the equation to allow the user to adjust theproportionality of transparency to pixel value. The final calculation ofalpha is: alpha=(source pixel/(maximum pixel value+1))^(alpha exponent)When the alpha exponent is 1.0, the relationship is linear. With valuesgreater than 1.0, the effect is to make the intermediate pixel valuesmore opaque toward the more grey. In any case, a 0 pixel value istransparent and a maximum pixel value is opaque.

The image data matrix resulting from the above process will be aprojected 2D view of a 3D image regenerated from the digitizedtomography view level image stack.

Features of the invention include:

-   1. Use of digitized tomography views to generate a 3D projection.-   2. Use of ray tracing to produce a projected view of a digitized    tomography view by using the intercept of a line and a plane, a line    from the view point through the projection plane cell, and a plane    representing a digitized tomography view level plane.-   3. Use of digitized tomography between pixel calculation to extract    value from projected view.-   4. Use of alpha blending to combine projected views.-   5. Using a function of pixel intensity to compute the alpha for    alpha blending.-   6. Use of an exponent to adjust the weighting of alpha in alpha    blending method:    -   a. increasing exponent increases opacity of ever lower pixel        values;    -   b. exponent of 1.0 same as no exponent;    -   c. can be used to adjust visible portion of image in projected        3D image.-   7. Increased computational efficiency by level shifting of view    point and projection point rather than shifting level of view level    plane:    -   a. shifting plane is computationally significant;    -   b. shifting projection point is combined with needed coordinate        transformation;    -   c. shifting view point is simple subtraction.

Conceptually, digitized tomography represents “analysis” in that itproduces the visual equivalent of a succession of “slices” of theexamined object. The 3D projection method added to digitized tomographyrepresents “synthesis” in that it produces the visual equivalent ofputting the “slices” back together. Therefore, the term digitizedtomosynthesis is appropriate. Digitized tomosynthesis has as its primaryvalue the ability to reduce visual complexity and highlight varioussubtle features. Similarly, the projection method as described can beused to project a reduced complexity 3D image by selecting only alimited range of view levels. Given the right sample and the right setupand configuration, the final result could offer the viewer theopportunity to see something that might not otherwise be obvious.

Math Model

The film mode simulation geometry works because the object plane isparallel to the film plane. Referring to FIG. 4, the geometry is suchthat the geometric magnification of projected shadows of the features ofthe object plane 18 is uniform. Thus the appropriate superposition ofthe film images places the projected shadows of the features of the viewimage plane in coincidence.

Referring to FIG. 5, in the alternate geometry, the object plane 18 isnot parallel to the film plane and thus the geometric magnificationvaries over the film plane. Since the film is dimensionally fixed, it istherefore not possible to superposition the film so that the projectedshadows of the features of the object plane are in coincidence. In orderto form an image by combining alternate geometry projections, one mustaccommodate the variation of magnification.

If one were to specify the geometry of image formation sufficiently, onecould mathematically trace a ray from the source 20, through the objectspace voxel to the image plane and compute to coordinate of the shadowof the voxel on the image plane. If done for each object rotation, theimage data could be extracted and combined to form the estimate of thegiven object space voxel. Since the ray tracing process includes thegeometric magnification effects, both the current and alternategeometries are applicable. As are many other geometries. The problem isin the specification of the geometry with sufficient accuracy andprecision to achieve acceptable results.

Limits to Resolution

There are three primary factors limiting resolution. The firstlimitation is the spot size of the x-ray source 20. Each edge of theprojected image has a penumbra caused by that spot size. Its dimensionis the spot size times the distance between the object feature and theimage plane divided by the distance between source 20 and the imageplane. Since the usual spot size is in the order of a millimeter, thesource distance in the order of a meter, and the feature distance is inthe order of a few tens of millimeters, the penumbra is in the order ofa few tens of thousandths of a millimeter or a few thousandths of aninch. The alternate geometry may have larger penumbra due to the largerdistances between the object feature and the image plane. Using asmaller spot size and/or longer source distances can reduce this effect.

The second limitation is the pixel size of the image plate. Since datafrom eight or more images is used, its possible to obtain resolutionsbetter than the pixel size of the image plate. Use of multiple imagesalso tends to cancel a portion of the penumbra effect.

The third limitation is the specification of the geometry. Errors in thespecification of the position of the source 20 have an effect similar inmagnitude as the penumbra effect. Those errors are attenuated by theratio of the source distance to feature distance. Errors in theidentification of the center of rotation has a one to one impact. Errorsin specification of the various angles are dependent upon the distanceover which the error is communicated. Acceptable errors for the currentimage plate are on the order of a few hundredths of a degree.

With such tight tolerance of angles and dimensions, the geometry must beeither very carefully measured or calibrated by using the system andspecially constructed objects.

Specifying the Geometry

The assumption is that there is a fixed relationship between the source20 and image plate, that the object is between the source 20 and theimage plate, and that it is the object that is rotated about a fixedaxis between each exposure cycle. If a means to measure the preciselocation and orientation of the source 20, the object, and image platewere to be developed, then it would be possible to remove thisrestriction. However, since the goal is to achieve a working systemwithin the bounds of currently available technology, the aboverestrictions are necessary. Relatively simple input transformations canbe used, if required, to adapt to precise and accurate location andorientation measurements

The Coordinate System

Looking at the system as it is presented in FIGS. 11 and 12, the X axisis left to right, the Y axis is vertical, and the Z axis is in and out.The origin of the system is the outer left most pixel on the imageplate. One coordinate increment is equal to one pixel. This makes theimage plate coincident with the X, Z plane. The X, Z plane in the newcoordinate system is the same as the X, Y plane of the original system.The change is necessitated by the three dimensional (3D) math requiredby the alternate geometries. The original geometry only required twodimensional (2D) math.

Dimensions to Specify

-   -   1. The source to image plate distance (source distance)    -   2. The image plate dimensions (width in Z axis, height in X        axis)    -   3. The image plate pixel size (from vendor)    -   4. The center of rotation to image plate distance (object        distance)    -   5. The object radius of interest (from operator)

The source to image plate distance is not critical. Errors arediminished by the ratio of the feature to image plate distance to thesource to image plate distance.

The image plate dimensions and pixel size are given by the manufacturer.

The center of rotation to image plate distance is also not critical. Itis used to compute the object plane level offset so that one can specifythe view level as distance from the base of the object.

The object radius of interest is perhaps the least critical. It simplysets the computed image size. It is necessary that it is sufficient toinclude the features of interest but not so large as to unnecessarilyextent the computation time.

Angles to Specify

-   -   1. The source to image plate plane angle in the YZ and XZ planes        (source angle)    -   2. The object axis of rotation to image plate plane angle in the        YZ and XZ planes (object angle)    -   3. Rotation angle of object for each image (rotation angle)

The source to image plate plane angles and the object axis of rotationto image plate plane angles are manufactured and/or measured as closelyas reasonably possible. Small errors can be measured using the systemand specially constructed objects and corrections can be computed. Ifthe manufacture or measurement is precise and accurate enough,corrections won't be required, in which case the ray tracing method isnot limited to any particular geometries.

It is apparent from our current system that the manufacture of the meansto set the rotation angle of the object is sufficiently precise andaccurate for the existing image plate.

Calibrating the Geometry

The existing system uses a two step geometry calibration. Referring toFIG. 6, the object 22 is a round marker disk 24 with a central pin 26made of a material suitably opaque to x-rays.

Referring to FIG. 7, in the first step, the disk 24 is retained on theobject turntable 28 via insertion of the bottom of the pin 26 throughthe disk 24 and into a hole in the center of the object turntable 28. Animage is taken and is used to discover the coordinates of the center ofrotation. In the alternate geometry, it can be used in the same manner

Referring to FIG. 8, in the second step, a plastic post 30 carries themarker disk and pin assembly at its top end. The plastic post 30 has asimilar pin at its bottom end, inserted into the central hole of theobject turntable 28, and a hole at its upper end to receive the bottomof the pin 26 through the round disk 24. A second image is taken and isused to discover the misalignment coordinate. The central voxel shadowcoordinate of the marker shadows is measured by the graphicalinteractive method, described above.

The center coordinate, the misalignment coordinate, the length of theplastic post 30, and the pixel size can be used to compute correctionsfor the source and object angles. This calibration technique is thesubject of my U.S. patent application Ser. No. 09/709,586, filed Nov.13, 2000, entitled “Tomosynthesis System and Registration Method,”hereby incorporated herein by reference. However, it was not thenextended to the alternate geometry. FIGS. 9 and 10 indicate how themarker disk 24 and misalignment post might be placed on an angled objectturntable 32 to calibrate the alternate geometry. The calibrationprocedure comprises the following steps

-   1. Compute the coordinates of an ordered set of object voxel    coordinates that form a foil-   2. Use that coordinate set to compute the corresponding image shadow    coordinates-   3. Extract and combine image shadow coordinates to form view foil    image pixels-   4. Assign view foil image pixels to corresponding pixels in an image    data array for presentation

A procedure similar to the one used with the film mode simulationgeometry could be used to compute and apply the corrections to thealternate geometry.

The alternate geometry may need one additional calibration step—thedistance between the image plate and the center of rotation (objectdistance). Errors in that dimension affect the accuracy of measurementof feature size. Simple measurement should be enough because the effectof that error is the size of error divided by the source distance.However, one way to calibrate that dimension is to place two marker disk24 and pin 26 assemblies a known distance apart on the object turntable32, centered on the center of the turntable 32 and oriented so the lineconnecting the centers of the disks is parallel to the Z axis. Thefollowing relationship would be true:

Solving for object distance can give a calibrated estimate of therequired dimension.

Ray Tracing Process

More particularly, the ray tracing process comprises the followingsteps:

-   1. Compute alignment calibration angles    -   a. Acquire centroid of center calibration marker shadow in image        plane        -   i. Coordinate is X1a,Z1a by method in Appendix C        -   ii. Note: the image plane is parallel with the XZ plane        -   iii. Note: the origin of the image plane is at the            coordinate system origin    -   b. Acquire centroid of misalignment calibration marker shadow in        image plane        -   i. Coordinate is X1b,Z1b by method in Appendix C        -   ii. Note: if Z1 a equals Z1b there is no misalignment    -   c. Compute misalignment angle in XZ plane        -   i. Misalignment angle equals arctan((X1b−X1a)/(z1b−z1a))−90°    -   d. Compute angle in XY plane        -   i. Convert height of misalignment calibration post to pixel            units        -   ii. Source angle equals arctan((misalignment post            height)/(X1b−X1a))        -   iii. For current geometry, angle is source angle            -   (1) Object angle is 0°            -   (2) Presumes accurate construction            -   (3) Implies object plane is parallel to image plane            -   (4) Construction and calibration procedure seems                adequate        -   iv. For alternate geometry, angle is object angle            -   (1) Source angle is 90°            -   (2) Presumes accurate alignment of image plane to source                path            -   (3) Implies optical axis is normal to image plane            -   (4) May need to develop a calibration procedure for                source angle-   2. Compute corrected 3D coordinates of source    -   a. Transform source distance to pixel units    -   b. Xs=X1a    -   c. Ys=(source distance)×sin(source angle)    -   d. Zs=Z1a-   3. Compute equation of image plate plane    -   a. Vector 1=(0,0,0),(0,(image width−1),0)        -   i. Vx1=0        -   ii. Vy1=image width−1        -   iii. Vz1=0        -   iv. M1=sqrt(Vx1^2+Vy1^2+Vz1^2)    -   b. Vector 2=(0,0.0),(0,(image height−1),0)        -   i. Vx2=0        -   ii. Vy2=image height−1        -   iii. Vz2=0        -   iv. MZ=sqrt(Vx2^2+Vy2^2+Vz2^2)    -   c. Normal Vector=the normalized cross product of Vector 1 and        Vector 2        -   i. Non-normalized vector            -   (1) NVx=(Vy1×Vx2)−(Vx1×Vy2)            -   (2) NVy=(Vz1×Vx2)−(Vx1×Vz2)            -   (3) NVz=(Vx1×Vy2)−(Vy1×Vx2)            -   (4) NM=sqrt(NVx^2+NVy^2+NVz^2)        -   ii. Normalized Vector            -   (1) NVx=NVx/NM            -   (2) NVy=NVy/NM            -   (3) NVz=NVz/NM    -   d. Plane coefficient        -   i. D=−NVx×Vx1−Nvy×Vy1−NVz×Vz1-   4. For each object voxel in desired digitized tomography view foil    -   a. Compute 3D coordinates of object voxel        -   i. Transform engineering units coordinate to pixel unit            coordinate        -   ii. Two dimensional coordinate transformations in XZ plane            -   (1) Rotate about axis of rotation to misalignment angle            -   (2) Rotate about axis of rotation to incremental angle        -   iii. For alternate geometry            -   (1) Two dimensional coordinate transformation in XY                plane                -   (a) Rotate about center of rotation to object angle                -   (b) Translate to the center of rotation    -   b. Use results of 1 and 4a to create a 3D line        -   i. (Xs,Ys,Zs),(Xo,Yo,Zo)    -   c. Use results of 3 and 4b to compute intercept coordinates        -   i. Dx=Xo−Xs;        -   ii. Dy=Yo−Ys;        -   iii. Dz=Zo−Zs;        -   iv. mu=−(D+NVx×Xd+Nvy×Yd+NVz*Zd)/(NVx×Dx+Nvy×Dy+NVz×Dz)        -   v. Xi=Xs+mu×Dx        -   vi. Yi=Ys+mu×Dy        -   vii. Zi=Zs+mu×Dz    -   d. For each image        -   i. Use results of 4c to extract projected ray pixel value            -   (1) Method in Appendix B    -   e. Combine pixel values from 4d by method in Appendix A    -   f. Store result of 4e in destination digitized tomography view        data matrix        Appendix A—Variable Absorption Digital Image Combination        Review of Absorption Digital Image Combination:

The process of producing a computed tomographic view from multipledigitized x-ray images requires the combination of the registered pixelsfrom each of the digitized images. The usual method used is a simplesummation or averaging of the pixels to be combined. This method resultsin an image containing many confusing artifacts and visual blurring.

The original absorption model for combining eight images is as follows:

-   -   The minimum pixel value was assigned a transmission fraction of        0.5    -   The maximum pixel value was assigned a transmission fraction of        1.0    -   All other pixel values were assigned a transmission fraction        proportionately between 0.5 and 1.0    -   For each pixel location in the computed view    -   Start with a maximum pixel value    -   Successively multiply initial pixel value by the transmission        fraction of that same pixel location in each of the eight images

This can be accomplished using integer arithmetic as follows:

-   -   Combined pixel=maximum pixel value    -   For each of eight images in succession        -   image pixel=image pixel+maximum pixel value        -   image pixel=image pixel DIV 2        -   image pixel=image pixel−maximum pixel value        -   image pixel=image pixel*combined pixel        -   image pixel=image pixel DIV maximum pixel value        -   combined pixel=combined pixel+image pixel    -   Note: While floating point arithmetic could be used, integer        arithmetic is much faster on most computers    -   Note: While eight images were used in the present case, the        method could be adapted to other numbers of images        Variable Absorption Digital Image Combination:

While the fixed absorption model gives generally good results, othervalues of absorption can give improved results for some exams. Forexample, exams of very dense objects can benefit from higher absorptionswhile exams with low density areas of interest can benefit from lowerabsorptions.

The fixed absorption model has been modified to allow variableabsorption using integer arithmetic as follows:

-   -   Combined pixel=maximum pixel value    -   For each of eight images in succession    -   image pixel=image pixel+maximum pixel value    -   image pixel=image pixel DIV 2    -   image pixel=image pixel−maximum pixel value    -   image pixel=image pixel*combined pixel    -   image pixel=image pixel DIV maximum pixel value    -   image pixel=image pixel*scaled absorption fraction    -   image pixel=image pixel DIV scaling factor    -   Combined pixel=combined pixel+image pixel

If the absorption fraction is one, the results are as before. Given anabsorption fraction greater than one, the effect is to increase theabsorption. An absorption fraction less than one results in a decreaseof absorption effect. Visual results are as expected.

Appendix B—High Resolution Digital Image Manipulation

Review of between pixel value calculation:

A digitized image is a rectilinear matrix of numbers. Each numberrepresents the integral intensity of a small area of the originalobject, picture, or transparency. In order to transform the matrix to animage the eye can use, hardware must be provided such that each numbercan be presented as a small area of light or dark on a video display orprinting device.

When the eye looks at a small area of light, a perception of brightnessresults from an integration of area and intensity. That is a small areaof bright will look as bright as a slightly larger area of slightly lessbrightness. This understanding lead to the following method ofdetermining the value of a pixel whose coordinate fell between pixels.

Given:

A two dimensional matrix of pixel values for a source image

A fractional coordinate of a pixel value to determinecoordinate=x.fx,y.fy,

-   -   Where:        -   x,y=the integral coordinate        -   fx=the fractional x coordinate        -   fy=the fractional y coordinate

There are four nearest pixels in the source image

-   -   Upper left=x,y    -   Upper right=x+l,y    -   Lower left=x,y+1    -   Lower right=x+l,y+1

There are four areas joining at x.fx,y.fy which are proportional to

-   -   Upper left=fx*fy*aspect ratio    -   Upper right=(1−fx)*fy*aspect ratio    -   Lower left=fx*(1−fy*aspect ratio)    -   Lower right=(1−fx)*(1−fy*aspect ratio)

Where:

-   -   aspect ratio=y pixel resolution/x pixel resolution

The eye will weight the four nearest pixels proportional to the opposingareas

Between pixel value=(Upper left value*lower right area+Upper rightvalue*lower left area+Lower left value*upper right area+Lower rightvalue*upper left area)/aspect ratio

Applications of between pixel value calculation:

The between pixel value calculation is currently used for twoapplications in the DT33 and DT36 programs:

-   -   Rotation and translation of digitized images    -   Forming the measure feature profile line

Tests have been performed yielding very good results for two moreapplications:

-   -   Digitized image magnification    -   Extended high resolution reference mark find

Several additional applications are possible:

-   -   Variable image area selection    -   Calculation of arbitrary z-axis view    -   Correction of digital image distortion    -   Magnification    -   Aspect ratio    -   Non-linear distortion

Review of rotation and translation of digitized images:

The digitized images must be rotated and translated for two reasons:

-   -   1 The acquired image must be rotated and translated to a        standard reference rotation and position    -   2 The set of images for an exam must be rotated and translated        so that the image combination will yield the computerized        tombgraphic view at the desired level

Given:

-   -   a two dimensional matrix of pixel values for a source image,    -   a two dimensional matrix of pixel value locations for a        destination image,    -   a coordinate of the image center of rotation,    -   a reference coordinate for the image center of rotation,    -   an angle to rotate about the image center of rotation,    -   an x displacement of the image center of rotation,    -   a y displacement of the image center of rotation,    -   transform each integral coordinate in the destination image into        a fractional coordinate in the source image:        -   x.fx=(x destination−source center x−x            displacement)*COS(rotation angle)+(y destination−source            center y+y displacement)*SIN(rotation angle)−(2.0*source            center x−reference x)        -   y.fy=(x destination−source center x−x            displacement)*SIN(rotation angle)+(y destination−source            center y+y displacement)*COS(rotation angle)+(2.0*source            center y−reference y)

Compute the between pixel value in the source image at x.fx,y.fy andplace it in the destination pixel location x destination,y destination.

Review of forming the measure feature profile line:

The operator selects two points on the image which are the ends of aline passing through the feature of interest.

Given the coordinates of the two points xl,yl and x2, y2,

Compute the slope of the given line:

-   -   IF x1<>x2 THEN slope=(y2−y1)/ (x2−x1)    -   ELSE slope=sign of (y2−y1)1E20 or similar large number

Compute the intercept of the given line:

-   -   intercept=y1−slope*x1

For each one pixel distance between the two points on the line computethe fractional coordinate:

-   -   Compute:        -   x increment=SQRT(1/(1+SQR(slope)))        -   y increment=slope*x increment        -   x.fx=x1+x increment*number of pixels from x1, y1        -   y.fy=y1+y increment*number of pixels from x1,y1

For each fractional coordinate compute the between pixel value to be thecorresponding feature profile line pixel value.

Digitized image magnification:

Previous applications were based upon the fractional coordinates beingone image pixel apart. If the fractional coordinates are spaced closerthan one pixel apart the destination image can be effectively magnified.That is if they are ½ pixel apart the effective magnification would be2×; ¼ pixel apart would be 4× and so on.

When such a magnification is done the result is a smooth magnificationof the digitized image-that appears much like the image produced by amagnifying glass. The customary approach to the magnification of adigital image is the simple square replication of each pixel. That isfor a 2× magnification each pixel is replicated into a 2 by 2 square.The resulting image appears magnified but has a square tile effect whichreduces the eyes ability to see the shape of the magnified features.

Magnification of digital images is becoming more important as theresolution of image digitizers and digital display devices improves.Current resolution is such that it is difficult to see features that area few pixels in size. This trend is expected to continue.

Extended high resolution reference mark find:

The current method of locating the center of the registration referencemarks is as follows:

The operator points to the inside of the reference mark.

The computer scans on the x and y axis until it finds a maximum in theslope of the image intensity higher than the noise peaks.

The fractional coordinates of the maximum slope is determined.

The position of the maximum is determined by using the maximum slopealong with the preceding and following slope values to form a secondorder curve.

The position of the maxima of the second order curve is assumed to bethe position of the edge.

The trial center of the mark is then computed as the mean x,ycoordinate.

The computer then takes the trial center and similarly computes a secondtrial center.

The compute then takes the second trial center and computes the finalcenter coordinate.

Since the reference marks are ellipsoidal:

The center coordinate of both reference marks are determined.

The errors in angle and position are computed.

The rotation and position of the image is adjusted by the rotation andtranslation procedure.

The computer then determines the new center coordinates of the referencemarks.

The sequence is repeated until the errors in angle and position arereduced below a small amount.

An extended method can be as follows:

The operator points to the inside of the reference mark.

The computer scans on the x and y axis until it finds a maximum in theslope of the image intensity higher than the noise peaks.

The fractional coordinates of the maximum slope is determined.

The position of the maximum is determined by using the maximum slopealong with the preceding and following slope values to form a secondorder curve.

The position of the maxima of the second order curve is assumed to bethe position of the edge.

The trial center of the mark is then computed as the mean x,ycoordinate.

The computer then takes the trial center and similarly computes a secondtrial center.

The compute then takes the second trial center and computes the finalcenter coordinate.

Do the above for both reference marks.

Compute the angle and positional error.

Use the same method used to create a feature profile line to produce anx,y scan across the reference marks whose angle is adjusted for theangular error found above.

In a similar manner find the fractional coordinates of the maximum slopein the pixel intensity.

Compute the center coordinates for both reference marks as the meancoordinate of each.

Compute a revised angle and position error.

Repeat the second process until the computed angle and position errorchanges less than some small amount.

The advantage of this extended method is that the true error in angleand position of the image is determined without the time consumingmultiple rotation and translation adjustment. This becomes varyimportant as the resolution of the digitized image improves. Currentimages are 3.1 meg pixels with a possible 12.6 meg pixels soon. Thehigher resolution image will take too much time and memory to allow suchadjustment.

Variable image area selection:

The current method of image area selection is as follows:

The image is adjusted to standard position and rotation.

The operator selects a packing fraction:

-   -   1 to 1    -   4 to 1    -   9 to 1    -   16 to 1.

The central area of the source image is packed into the destinationimage area.

This method has several problems:

The time and memory required to adjust the digitized image to standardposition and rotation.

The inflexibility of the integral pixel packing approach not allowingoptimum sizing of the area of interest after digitization.

Demagnification of the digitized image should be possible in a mannervary similar to the magnification discussed above. Compute thefractional coordinates so that they are more than a source pixel apart.Any arbitrary demagnification should be able to be achieved.

Rotation and translation could be accomplished at the same time by usingthe extended reference mark find method coupled with the coordinatetransformation procedure above. Thus the excessive time and memoryrequired for the current method could be eliminated along with allowingan optimum selection of the area of interest.

Calculation of Arbitrary Z-Axis View:

The current method of computation of a z-axis view is limited to any xor y line in the currently formed computed tomography view series. Aspecific x or y line is taken from each stored computed tomographic viewand displayed side by side. This is not quite satisfactory because:

-   -   The entire series must be created first.    -   Too much time is taken read all of the data from the disk.    -   The feature of interest may not be along the x or y axis.    -   The feature of interest may in fact be curved.        An alternative can be:

The operator selects an arbitrary path on one view.

A list of fractional coordinates one pixel apart are generated from thatpath or possibly magnified or demagnified as required.

Transform that list so that the appropriate rotations and translationsare included.

Acquire pixel data along that transformed path by the above betweenpixel method.

Combine the pixels by the currently used absorption method.

Present the computed paths side by side.

Correction of digital image distortion:

If an array of fractional coordinates can be computed to correct for thedistortion in a digital image, the above between pixel calculation canbe used to correct that distortion. Thus linear distortion such asmagnification or aspect ratio can be corrected with ease. Non-lineardistortion such as that resulting from a variable speed scanning cameracould be corrected if that non-linearity were adequately known.

Appendix C—High Resolution Registration of Digital Images

Introduction

The coordinate system to be used in this disclosure is as follows:

-   -   The film is placed at the image plane.    -   The image plane is the x, y plane.    -   The image plane is tilted at some angle with respect to a base        plane.    -   The y axis is parallel to be base plane.    -   The image is formed by x-rays projecting a shadow of an object        at or near the image plane.    -   Multiple images are formed by rotating the object about an axis        of rotation    -   The z axis is parallel to the axis of rotation.    -   The intersect of the axis of rotation and the image plane is the        center of rotation.    -   The center of rotation is the desired origin of the coordinate        system.

Previous disclosures described a method of using reference marks totransfer the registration of the object to the computer. Several methodsof locating and using the centers of those markers were discussed.

As the image resolution is pushed to higher levels, registration becomesa more demanding problem to be solved. For digitized tomography brandcomputerized tomography to work at its best, the multiple digitizedimages must be in registration to within a pixel.

After digitization of the x-ray images, the two registration markscenters are located, the apparent center of rotation is computed basedupon a previously established relationship between the registrationsmarks and the center of rotation. The image is then moved to a standardposition and rotation:

-   -   The line joining the centers of the registration marks is        parallel to the y axis.    -   The center of rotation is positioned at the center of the        digital image.

Since the acquired images must be rotated to form the final CT view, thequality of the result is vary sensitive to variations in registrationalong the y axis. This is especially true as the area of interest ismore distant from the center of rotation.

When working with pixels as small as a few mills, the projected image ofthe reference marks becomes an ellipse rather than a circle. Center findroutines based upon the assumption of a circle fail to locate the centerof the mark to within the required fraction of a pixel. In fact,variations in the position of the final image of more than five pixels,have been observed.

Methods have been developed that markedly reduces this problem. Reviewof previously disclosed center find methods:

Step Description

-   1 Locate the coordinates of all reference mark edge pixels-   2 Compute mean x and y coordinate-   3 Use computed x, y coordinate as center of reference mark    Method One:    Step Description-   1 Use mouse to point to inside of reference mark.-   2 Use pointed to pixel as starting point.-   3 Scan along the x axis to find maxima or minima of differential of    pixel level.-   4 Compute x coordinate of left and right edge of reference mark as    position of maxima or minima of differential.-   5 Compute mean x coordinate.-   6 Scan along the y axis to find maxima or minima of differential of    pixel level.-   7 Compute y coordinate of top and bottom edge of reference mark as    position of maxima or minima of differential.-   8 Compute mean y coordinate.-   9 Use computed mean x,y coordinate as center of reference mark    Method two—in use:    Method Two:    Step Description-   1 Use mouse to point to inside of reference mark.-   2 Use mouse to point to outside of reference mark.-   3 Use pointed to inside pixel as starting point.-   4 Use a threshold pixel level based upon the pointed to inside and    outside pixel levels.-   5 Scan along the x axis until the first pixel equal to or less than    the threshold level is found.-   6 Use the x coordinate of the found pixel as the right and left edge    x coordinate.-   7 Compute the mean x coordinate.-   8 Scan along the y axis until the first pixel equal to or less than    the threshold level is found.-   9 Use the y coordinate of the found pixel as the top and bottom edge    y coordinate.-   10 Compute the mean y coordinate.-   11 Use computed mean x,y coordinate as the center of reference mark    Improvements being tested.    Method Two Improvement:

Instead of using the coordinate of the pixel equal to or less than thethreshold level, use the interpolated factional coordinate:

-   -   Left coordinate=x coordinate+((pixel x+1        level)−threshold)/((pixel x+1 level)−(pixel x level))    -   Right coordinate=x coordinate−((pixel x        level)−threshold)/((pixel x level)−(pixel x−1 level))    -   Top coordinate=y coordinate+((pixel y+1        level)−threshold)/((pixel y+1 level)−(pixel y level))    -   Bottom coordinate=y coordinate−((pixel y        level)−threshold)/((pixel y level−(pixel y−1 level))

This improvement changes method two resolution to much better than itsprevious one half pixel.

Improved estimation of center of reference mark:

Using either method one or two locate the approximate center of thereference mark. Use the first found value of the center as the startingpoint of a second approximation of the center.

This improvement changes the reproducibility of the found center tobetter than one tenth of a pixel. Positioning by method of successiveapproximation:

Step Description

-   0 Use one of the previously disclosed methods of finding the    centers.-   1 Compute the angle error and center of rotation of the image.-   3 Adjust the image to an approximate standard position and rotation.-   4 Find the new centers of the reference marks using the same find    method.-   5 Compute the new angle error and center of rotation.-   6 If the new angle error and center of rotation are too large go to    step 3.

This improvement changes the reproducibility of the rotation andposition of the image to better than five hundredths of a degree and onepixel.

This method can be performed using the current software. However, itrequires repeated operator interaction and a substantial amount of timefor the repeated adjustment to standard position and rotation.

Step Description

-   1 Find the centers of the reference marks by either improved method.-   2 Use the right, left, top, and bottom coordinates to determine the    size of the reference mark.-   3 Use the same method used in the measure image feature to create    the pixel pattern across the reference marks as follows:    -   a: The diameter plus a small amount through the center parallel        to the line between centers.    -   b: The diameter plus a small amount through the center        perpendicular to the line between centers.-   4 Determine the maximum and minimum pixel level in the pixel    patterns.-   5 Compute the edge coordinates of the half way level in each    pattern.-   6 Use the mean coordinates as the center coordinates for step 3    until the difference between the new center and the previously found    center is less than a to be specified fraction of a pixel.-   7 Use the final found centers as a basis for adjusting the digital    image to standard position and rotation.

This method can have the advantage of requiring only one sequence ofoperator actions and only one adjustment to standard position androtation.

Review of the measure image pixel pattern:

Step Description

-   1 Determine two different position coordinates on the image.-   2 Compute the equation of the line between the two coordinates.-   3 Compute the image coordinates of one pixel increments along that    line.-   4 Determine the pixel value at each pixel pattern coordinate:    -   a: The pixel pattern coordinates may be fractional coordinates.    -   b: The image pixels are at integral coordinates.    -   c: The pixel value will be determined by the same method used by        the rotate—translate routine—an opposing area weighted mean of        the nearest four pixel values.        Appendix D—Axial Misalignment Calibration for Non-Film X-Ray        Device Images

The digitized tomography process requires very precise imageregistration for it to produce high quality computed images. Therequirement is such that every pixel in every source image must be knownto the process to within a fraction of a pixel. A method has beendeveloped that provides for such high resolution image registration offilm x-ray images. It consisted of a fixture, a setup process, andcalibration procedure that assured a very specific geometry of exposureand allowed each digitized film image to be transformed into a standardregistration and format. It was specifically because the x-ray sourcewas position immediately overhead the center of rotation of the fixtureand the fixture was tilted to the desired angle that axial misalignmentof the x-ray images was not encountered.

The current digitized tomography calculations assume the path betweenthe center of the x-ray source and the center of rotation of thenon-film fixture's turntable 28 falls on a plane that is perpendicularto the non-film image plate's surface, parallel to the long axis (the“Y” axis) of the non-film image plate, and passes through the center ofrotation of the fixture's turntable 28. Since the current fixture isplaced horizontally and the x-ray source is moved to the side and angledtoward the center of rotation, an axial misalignment is possible withoutusing special alignment procedures and tools. Even a small misalignmentcan cause serous degradation of the resultant computed digitizedtomography views.

Currently, the center of rotation is calibrated using an image of acircular marker placed at the center of rotation in a manner identicalto the method used to identify the center of rotation of film images.Because all images can be taken with an identical registration, thethree reference mark technique used with film images is not required.However, a second reference mark must be used to measure the axialmisalignment of the x-ray optical axis. That marker may be of the samesize and material as the center reverence mark but must be displaceabove the center of rotation by some distance. Since the center positionof the marker can be determined to approximately ½ of a pixel width,that distance need only be approximately ½ the maximum possible imageradius. The current image plate is 2304 pixels wide with a 200 pixel“gutter” along one edge. That makes the effective radius(2304-200)/2=1145 pixels at 0.005 inches per pixel, that would be 5.724inches. One half of that would be 2.86 inches.

Using the center coordinates of the two reference markers, the axialmisalignment angle can be determined and a correction to the digitizedtomography calculations can be applied to the “X” and “Y” axisdisplacement values. The process is as follows:

Let x1,y1 be the pixel coordinate of the center of the center referencemark image. Let x2,y2 be the pixel coordinate of the center of theelevated center reference mark image. Let the x-ray optical axis be moreor less parallel to the “Y” axis in the positive direction. The distancebetween centers of the reference marks can be computed by taking thesquare root of the sum of the squares of the differences between thefirst and second coordinates:distance=sqrt((x1−x2)*(x1−x2)+(y1−y2)″(y1−y2)). That distance can bedivided into the differences between the coordinates to compute the sineand cosine of the misalignment angle: sine of misalignmentangle=(x2−x1)/distance and cosine of misalignmentangle=(y2−y1)/distance.

The displacement is the normally calculated displacement value and iscomputed as a function of the image level. The “X” axis displacementvalue can be computed as the displacement times the sine of themisalignment angle and the “Y” axis displacement value can be computedas the displacement times the cosine of the misalignment angle.

The adjusted displacements are applied to their respective axiscalculations as described elsewhere.

Appendix E—Basic 3D Math

This article is part of The Win95 Game Programmer's Encyclopedia

Basic 3D Math

Introduction

The Win95GPE article on basic 2D math lays a lot of the ground work forunderstanding the information in this article. This section containssome of the basic high-school math needed to understand the equationsand graphics primitives used in 3D computer graphics. 3D math is oftenjust a simple extension of the same techniques used for 2D math. An aidto understanding 3D math is to think of it as simply another programminglanguage.

Definition of a 3D Point

A point is similar to it's 2D counterpart, we simply add an extracomponent, Z, for the 3rd axis, as shown in FIG. 11. Points are nowrepresented with 3 numbers: <x,y,z>. This particular method ofrepresenting 3D space is the “left-handed” coordinate system, which Iuse throughout this entire article. In the left-handed system the x axisincreases going to the right, the y axis increases going up, and the zaxis increases going into the page/screen. The right-handed system isthe same but with the z-axis pointing in the opposite direction.

Distance Between Two 3D Points

The distance between two points <Ax,Ay,Az> and <Bx,By,Bz> can be foundby again using the pythagorus theorem:dx=Ax−Bxdy=Ay−Bydz=Az−Bzdistance=sqrt(dx*dx+dy*dy+dz*dz)Definition of a 3D Vector

Like it's 2D counterpart, a vector can be thought of in two ways: eithera point at <x,y,z> or a line going from the origin <0,0,0> to the point<x,y,z>.

3D Vector addition and subtraction is virtually identical to the 2Dcase. You can add a 3D vector <vx,vy,vz> to a 3D point <x,y,z> to getthe new point <x′,y′,z′> like so:x′=x+vxy′=y+vyz′=z+vz

Vectors themselves can be added by adding each of their components, orthey can be multiplied (scaled) by multiplying each component by someconstant k (where k< >0). Scaling a vector by 2 (say) will still causethe vector to point in the same direction, but it will now be twice aslong. Of course you can also divide the vector by k (where k< >0) to geta similar result.

To calculate the length of a vector we simply calculate the distancebetween the origin and the point at <x,y,z>:$\begin{matrix}{{length} = {{{< x},y,{z > {- {< 0}}},0,{0 >}}}} \\{{sqrt}\left( {{\left( {x - 0} \right)*\left( {x - 0} \right)} + {\left( {y - 0} \right)*\left( {y - 0} \right)} + {\left( {z - 0} \right)*\left( {z - 0} \right)}} \right)} \\{{sqrt}\left( {{x*x} + {y*y} + {z*z}} \right)}\end{matrix}$

Often in 3D computer graphics you need to convert a vector to a unitvector, i.e., a vector that points in the same direction but has alength of 1. This is done by simply dividing each component by thelength:Let <x,y,z> be our vector, length=sqrt(x*x+y*y+z*z)${{Unit}\quad{vector}} = {\frac{{< x},y,{z >}}{length} = {{\frac{x}{length},\frac{y}{length},\frac{z}{length}}}}$(where length=|<x,y,z>|)

Note that if the vector is already a unit vector then the length will be1, and the new values will be the same as the old.

Definition of a Line

As in 2D, we can represent a line by it's endpoints (P1 and P2) or bythe parametric equation:P=P1+k*(P2−P1)where k is some scalar value between 0 and 1.The Dot Product

The dot product between two vectors <Ax,Ay,Az> and <Bx,By,Bz> iscalculated like so:A*B=Ax*Bx+Ay*By+Az*Bz

If A and B are unit vectors then the dot product is the cosine of theangle between them, so the angle itself can be calculated by taking theinverse cosine of the dot product:theta=invcos(A*B)

Fortunately you'll often only need the cosine of the angle between 2vectors and not the angle itself, so the expensive step of calculatingthe inverse cosine can be skipped.

Definition of a Plane

A plane is an infinitely wide flat polygon-type surface. A plane can bedefined by a normal <nx,ny,nz> and a scalar value k. The normal can bethought of as representing the direction that the surface of the planeis facing, i.e., it's a directional vector at a right angle to theplane, as shown in FIG. 12. We can imagine the normal as describing aninfinitely long line stretching off to infinity in either direction, andthat a plane can slide up and down along this line. In the planeequation the value k specifies where exactly on this line the planesits.

The equation for the plane uses the dot product:<x,y,z>*<nx,ny,nz>=k

All points <x,y,z> that actually lie on the plane itself will satisfythis equation. This gives us a convenient method for determining whichside of a plane any given point <x,y,z> is: $\begin{matrix}{{< x},y,{z > {*{< {nx}}}},{ny},{{{nz} >}\quad = k}} & {{Point}\quad{is}\quad{in}\quad{plane}} \\{{< x},y,{z > {*{< {nx}}}},{ny},{{nz} > \quad > k}} & {{Point}\quad{is}\quad{on}\quad{one}\quad{side}\quad{of}\quad{plane}} \\{{< x},y,{z > {*{< {nx}}}},{ny},{{nz} > < k}} & {{Point}\quad{is}\quad{on}\quad{other}\quad{side}\quad{of}\quad{plane}}\end{matrix}$

The vector <nx,ny,nz> and scalar k are unique to every plane (a way ofcalculating them is shown below). These equations are helpful inperforming back-face culling. Substitute the view point into theequation, if the value comes out less than k then you know that you arefacing the “back” side of the polygon and thus don't need to draw it.

The Cross Product

If you have 2 vectors then the cross product will return a vector whichis perpendicular (i.e., at a right angle).to both of them. The crossproduct between two vectors <Ax,Ay,Az> and <Bx,By,Bz> is:A×B=<Ay*Bz−Az*By, Az*Bx−Ax*Bz, Ax*By−Ay*Bx>

Note that while the dot product returns a single scalar value, the crossproduct returns a vector.

Taking the cross-product between the x axis <1,0,0> and y axis <0,1,0>gives us: $\begin{matrix}{{< 1},0,{0 > x < 0},1,{0>= < {{0*0} - {0*1}}},{{0*1} - {1*0}},{{{1*1} - {0*0}} >}} \\{{= {< 0}},0,{1 >}}\end{matrix}$which is of course the z axis and is indeed perpendicular to bothvectors.Calculating a Plane from 3 Points

To calculate a plane from 3 given points we first calculate the normal.If we imagine the 3 points form three edges in the plane then we cantake two of the edges and calculate the cross-product between them. Theresulting directional vector will be the normal, and then we can plugany of the 3 known points into the plane equation to solve for k. Forpoints p1,p2 and p3 we get:normal=(p1−p2)×(p3−p2)k=normal*p1as shown in FIG. 13. Note that it is extremely important to keep trackof which direction your points are stored in. Take 3 points stored inclockwise direction in the x/y plane, as shown in FIG. 14. The normal tothe plane these 3 points define is: $\begin{matrix}{{normal} = {\left( {{p1} - {p2}} \right) \times \left( {{p3} - {p2}} \right)}} \\{= {\left( {0,{- 1},0} \right) \times \left( {1,{- 1},0} \right)}} \\{{= {< {{\left( {- 1} \right)*0} - {0*\left( {- 1} \right)}}}},{{0*1} - {0*0}},{{{0*\left( {- 1} \right)} - {\left( {- 1} \right)*1}} >}} \\{{= {< 0}},0,{1 >}}\end{matrix}$i.e., the z axis. If we were to store the points counter-clockwise thenormal calculated would be<0,0,−1>, which is still the z axis but in the“opposite” direction. It's important to keep track of these things sincewe often need plane equations to be correct in order to determine whichside of a polygon an object (such as the view point) is on.3D Rotation

In the Win95GPE article on basic 2D math we saw how a 2D point <x,y> canbe rotated around the origin <0,0>. In 3D this is the same as rotatingthe point around the z axis (the z value remains the same). We canslightly modify this basic equation to get rotation around all threeaxis: $\begin{matrix}\underset{\_}{{Rotation}\quad{about}\quad{the}\quad x\quad{{axis}:}} \\\begin{matrix}{x^{\prime} = x} \\{y^{\prime} = {\left( {\cos\quad é*y} \right) - \left( {\sin\quad é*z} \right)}} \\{z^{\prime} = {\left( {\sin\quad é*y} \right) + \left( {\cos\quad é*z} \right)}}\end{matrix} \\\underset{\_}{{Rotation}\quad{about}\quad{the}\quad y\quad{{axis}:}} \\\begin{matrix}{x^{\prime} = {\left( {\cos\quad é*x} \right) + \left( {\sin\quad é*z} \right)}} \\{y^{\prime} = y} \\{z^{\prime} = {{- \left( {\sin\quad é*x} \right)} + \left( {\cos\quad é*z} \right)}}\end{matrix} \\\underset{\_}{{Rotation}\quad{about}\quad{the}\quad z\quad{{axis}:}} \\\begin{matrix}{x^{\prime} = {\left( {\cos\quad é*x} \right) - \left( {\sin\quad é*y} \right)}} \\{y^{\prime} = {\left( {\sin\quad é*x} \right) + \left( {\cos\quad é*y} \right)}} \\{z^{\prime} = z}\end{matrix}\end{matrix}$

Rotating a point around an arbitrary axis is a tad complex, so I'lldiscuss how to do this in the article on matrices (since they're anideal tool for performing these types of rotations). It's fairly easy toexpand the steps in that article for “regular” 3D math.

Intersections

This section shows how to calculate intersections between variousobjects. This is particularly handy for things like collision detection.

Intersection Between a Line and a Plane

This occurs at the point which satisfies both the line and the planeequations.

 Line equation: p=org+u*dir  (1)Plane equation: p*normal−k=0.  (2)

Substituting (1) into (2) and rearranging we get:(org+u*dir)*normal−k=0i.e., u*dir*normal=k−org*normali.e., u=(k−org*normal)/(dir*normal)

If (d*normal)=0 then the line runs parallel to the plane and nointersection occurs. The exact point at which intersection does occurcan be found by plugging u back into the line equation in (1).

Intersection Between a Line and a Sphere

This occurs at the point which satisfies both the line and the sphereequations.Line equation: p=org+u*dir  (1)Sphere equation: |p−origin|=radius  (2)

Substituting (1) into (2) we get:|(org+u*dir)−origin|=radiusi.e., (org+u*dir−origin)^2=radius^2

Which can be rearranged into the following form:u^2*dir^2+u*2*dir*(org−origin)+(org−origin)^2

This is a quadratic equation in u which can thus be solved with thefollowing formula:u=−B+/−sqrt(B^2−4AC)2Awhere A=dir^2

-   -   B=2*dir*(org−origin)    -   C=(org−origin)^2

Note that dir^2 means dir*dir, i.e., the dot product of dir with itself.The same applies to org−origin in C. dir*(org−origin) in B is also a dotproduct.

To get the actual points of intersection you plug the 2 resulting uvalues into the line equation.

If A=0 then the line does not intersect the sphere. If sqrt(B^2−4AC)=0then the line is tangent to the surface of the sphere and there is onlyone point of intersection at u=−B/2A.

Intersection Between Three Planes

In order to find the point <x,y,z> of intersection between three planes(p*n1−k1=0, p*n2−k2=0 and p*n3−k3=0) we can use matrices (Appendix F; inparticular read the section on solving simultaneous equations). If weplace all three equations into a matrix then we can see that thefollowing rule will apply: ${{\begin{matrix}{{n1}.x} & {{n1}.y} & {{n1}.z} & {- {k1}} \\{{n2}.x} & {{n2}.y} & {{n2}.z} & {- {k2}} \\{{n3}.x} & {{n3}.y} & {{n3}.z} & {- {k3}} \\0 & 0 & 0 & 1\end{matrix}} \times {\begin{matrix}x \\y \\z \\1\end{matrix}}} = {\begin{matrix}0 \\0 \\0 \\1\end{matrix}}$

By rearranging this equation we can solve for <x,y,z>: ${\begin{matrix}x \\y \\z \\1\end{matrix}} = {{\begin{matrix}{{n1}.x} & {{n1}.y} & {{n1}.z} & {- {k1}} \\{{n2}.x} & {{n2}.y} & {{n2}.z} & {- {k2}} \\{{n3}.x} & {{n3}.y} & {{n3}.z} & {- {k3}} \\0 & 0 & 0 & 1\end{matrix}}^{- 1} \times {\begin{matrix}0 \\0 \\0 \\1\end{matrix}}}$

Thus if we calculate the inverse of the matrix then the 4th column willcontain the point of intersection. (If the matrix has no inverse then atleast two of the planes are parallel and there is no point ofintersection).

Note that because of the constant terms we could instead use a 3×3matrix for a more optimized solution: ${\begin{matrix}x \\y \\z\end{matrix}} = {{\begin{matrix}{{n1}.x} & {{n1}.y} & {{n1}.z} \\{{n2}.x} & {{n2}.y} & {{n2}.z} \\{{n3}.x} & {{n3}.y} & {{n3}.z}\end{matrix}} \times {\begin{matrix}{k1} \\{k2} \\{k3}\end{matrix}}}$Intersection Between Two Planes

The line (p=org+u*dir) of intersection between two planes (p*n1−k1=0 andp*n2−k2=0) is perpendicular to each plane's normal, i.e.:dir=n1×n2

What we can now do is assume there is a third plane perpendicular toboth planes, i.e., one which has a normal of dir. All that remains is tofind the k value for this plane and then find the point of intersectionbetween all three planes, as shown in the previous section. Forsimplicity we can assume this third k value to be 0 and use theresulting point as our org value: ${{\begin{matrix}{{org}.x} \\{{org}.y} \\{{org}.z} \\1\end{matrix}} = {{\begin{matrix}{{n1}.x} & {{n1}.y} & {{n1}.z} & 0 \\{{n2}.x} & {{n2}.y} & {{n2}.z} & 0 \\{{dir}.x} & {{dir}.y} & {{dir}.z} & 0 \\0 & 0 & 0 & 1\end{matrix}} \times {\begin{matrix}{k1} \\{k2} \\0 \\1\end{matrix}}}}\quad$

Copyright (c) 1997 Mark Feldman (pcgpe@geocities.com)—All RightsReserved

This article is part of The Win95 Game Programmer's Encycloopedia

Please retain this footer if you distribute this file

Appendix F—Matrices

This article is part of The Win95 Game Programmer's Encyclopedia

Matrices

Introduction

Matrices are extremely handy for writing fast 3D programs. They are justa 4×4 list of numbers, but they do have 2 very important properties:

-   -   1) They can be used to efficiently keep track of        transformations, i.e., actions which occur in a VR program such        as movement, rotation, zoom in/out etc.    -   2) A single matrix can represent an infinite number of these        transformations in any combination. Let's say the user in your        program walks forward, turns left, looks up, backs up a bit etc        . . . All you need to do is keep a copy of a master matrix in        memory and adjust it as the user does these things. At any point        you then use this one matrix to figure out where everything in        your virtual world should be drawn on the screen.

A transformation is simply a way of taking a set of points and modifyingthem in some way to get a new set of points. For example, if the usermoves 10 units forward in a certain direction then the net result is thesame as if all objects in the world moved 10 units in the oppositedirection.

In 3D computer graphics it is often convenient to assume that theview-point is always at the origin <0,0,0>. Not only does it make iteasier to render a scene but it often makes it a lot faster as well.Instead of having the user move around the world we can keep the user atthe origin and make the objects move around instead. From the userspoint of view it will look exactly the same. Matrices are a fast andconvenient tool to perform these transformations.

A Point in Space

As mentioned above, a point or vector can be represented in 3D spaceas<x,y,z>. When using matrix math it helps to represent it as<x,y,z,w>.That extra w there helps make things like movement easier to deal with.If we have a point in 3D space we can convert it to this new format bysetting w to 1. In this text I'll be representing all points in space asa column vector:${{\begin{matrix}x \\y \\z \\w\end{matrix}} < {- w}} = 1$Modifying the Position of a Point

Let's say we want to take any point <x,y,z,w> given and do something toit to get a new point. A row vector can be used to represent how we wantto change a point: $\begin{matrix}{{These}\quad{values}\quad{contain}\quad{the}\quad{info}\quad{on}} & \left. \longrightarrow\quad \right. & {\begin{matrix}A & B & C & D\end{matrix}} & \cdot & {\begin{matrix}x \\y \\z \\w\end{matrix}} & \longleftarrow & {{This}\quad{is}\quad{our}\quad{point}\quad{in}\quad 3D\quad{space}} \\{\quad{{how}\quad{we}\quad{want}\quad{to}\quad{change}\quad{the}\quad{point}}\quad} & \quad & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$

To figure out how the row vector |A B C D| changes a point, we canvisualize lying the point's column vector on it's side on top of it likethis:${\begin{matrix}x & y & z & w \\A & B & C & D\end{matrix}} = {\left( {x*A} \right) + \left( {y*B} \right) + \left( {z*C} \right) + \left( {w*D} \right)}$

Compare this to the article on basic 3D math (Appendix E) and you'll seethat we are in fact taking the dot product of the two vectors. What wedo above is multiply each top item by the item under it and add theresults up to get the answer.

Let's say we want to be able to take any point <x,y,z,w> and figure outthe coordinates for the point <x′,y′,z′,1> which is exactly 4 units tothe “right” of it, i.e., the point which is 4 units further along the xaxis. We can do this by using 4 row vectors. The first one willcalculate the new x point (x′), the next one the new y point (y′) and soon. First let's look at calculating x′.

We know that the new x point can be calculated like this: x′=x+4, so arow vector to calculate this would be:|1 0 0 4|i.e. when we multiply this out by a point <x,y,z,w> we'll get:(x*1)+(y*0)+(z*0)+(w*4)=x+4

We also know that y′=y, z′=z and w=1, so we can calculate the rowvectors for each of the other values, and stack them all on top of eachother: $\begin{matrix}x & {{row}\quad{vector}} & \rightarrow & \quad & 1 & 0 & 0 & 4 & \quad & {x^{\prime} = {{{1*x} + {0*y} + {0*z} + 4} = {x + 4}}} \\y & {{row}\quad{vector}} & \rightarrow & \quad & 0 & 1 & 0 & 0 & \quad & {x^{\prime} = {{{0*x} + {1*y} + {0*z} + 0} = y}} \\z & {{row}\quad{vector}} & \rightarrow & \quad & 0 & 0 & 1 & 0 & \quad & {x^{\prime} = {{{0*x} + {0*y} + {1*z} + 0} = z}} \\1 & {{row}\quad{vector}} & \rightarrow & \quad & 0 & 0 & 0 & 1 & \quad & {x^{\prime} = {{{0*x} + {0*y} + {0*z} + 1} = 1}}\end{matrix}\quad$

And VOILA! That's what a matrix is! To take a point <x,y,z,w> andcalculate the new point we just multiply the point by each of the rowvectors.

Here's a more generic representation of what we are doing:${\begin{matrix}x^{\prime} \\y^{\prime} \\z^{\prime} \\w^{\prime}\end{matrix}} = {{\begin{matrix}{M11} & {M12} & {M13} & {M14} \\{M21} & {M22} & {M23} & {M24} \\{M31} & {M32} & {M33} & {M34} \\{M41} & {M42} & {M43} & {M44}\end{matrix}}\quad{\begin{matrix}x \\y \\z \\w\end{matrix}}}$

In this case <x,y,z,w> is the point we are popping in, <x′,y′,z′,w′> isthe point we'll be getting out, and each number in the matrix isrepresented by Mij, where i=row number and j=column number.

Following the line of reasoning above, we can figure out how tocalculate the new point based on the given point and the matrix:$\begin{matrix}{x^{\prime} = {\left( {x*{M11}} \right) + \left( {y*{M12}} \right) + \left( {z*{M13}} \right) + {M14}}} \\{y^{\prime} = {\left( {x*{M21}} \right) + \left( {y*{M22}} \right) + \left( {z*{M23}} \right) + {M24}}} \\{z^{\prime} = {\left( {x*{M31}} \right) + \left( {y*{M32}} \right) + \left( {z*{M33}} \right) + {M34}}} \\{w^{\prime} = {\left( {x*{M41}} \right) + \left( {y*{M42}} \right) + \left( {z*{M43}} \right) + {M44}}}\end{matrix}$

In practice, we don't really need that last equation, we know that itwill always equal 1 so we don't have to calculate it.

Plugging a point into a matrix like this and getting a new point out iscalled “transforming” a point by a matrix, just keep that in mind so youknow what to call your procedure to do it!

A few of the more commonly used matrices follow.

Doing Nothing

You often need a matrix to pop out exactly the same point as was pluggedin. The matrix to do this is called the identity matrix: $\begin{matrix}1 & 0 & 0 & 0 & \quad & \quad & \quad & \quad & {x^{\prime} = x} \\0 & 1 & 0 & 0 & \quad & \quad & \quad & \quad & {y^{\prime} = y} \\0 & 0 & 1 & 0 & \quad & \quad & \quad & \quad & {z^{\prime} = z} \\0 & 0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$Translation

A translation is simply adding (or subtracting) a point to any givenpoint. Let's say you want to add TX to x, TY to y and TZ to z. Thematrix to do this is: $\begin{matrix}1 & 0 & 0 & {tx} & \quad & \quad & \quad & \quad & {x^{\prime} = {x + {tx}}} \\0 & 1 & 0 & {ty} & \quad & \quad & \quad & \quad & {y^{\prime} = {y + {ty}}} \\0 & 0 & 1 & {tz} & \quad & \quad & \quad & \quad & {z^{\prime} = {z + {tz}}} \\0 & 0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$Scaling

Sometimes we may need to scale a point, i.e., mutiply each axis by agiven number. This is handy for things like zoom effects.$\begin{matrix}{sx} & 0 & 0 & 0 & \quad & \quad & \quad & \quad & {x^{\prime} = {{sx}*x}} \\0 & {sy} & 0 & 0 & \quad & \quad & \quad & \quad & {y^{\prime} = {{sy}*y}} \\0 & 0 & {sz} & 0 & \quad & \quad & \quad & \quad & {z^{\prime} = {{sz}*{tz}}} \\0 & 0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$

Of course, if you only want to scale along the X axis (say) then yousimply set SX to the scale value and set SY and SZ to 1.

Basic Rotations

We can rotate things around the x axis, the y axis, or the z axis. Eachaxis forms a line in 3D space. We can also rotate things around anyarbitrary line. This is handy if we want to do effects like objectsrotating about their own axis (which will be discussed later).

The matrix for rotating around the x axis by angle é is: $\begin{matrix}1 & 0 & 0 & 0 & \quad & \quad & \quad & \quad & {x^{\prime} = x} \\0 & {\cos\quad é} & {{- \sin}\quad é} & 0 & \quad & \quad & \quad & \quad & {y^{\prime} = {{\left( {\cos\quad é} \right)*y} - {\left( {\sin\quad é} \right)*z}}} \\0 & {\sin\quad é} & {\cos\quad é} & 0 & \quad & \quad & \quad & \quad & {z^{\prime} = {{\left( {\sin\quad é} \right)*y} + {\left( {\cos\quad é} \right)*z}}} \\0 & 0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$

The matrix for rotating around the y axis by angle é is: $\begin{matrix}{\cos\quad é} & 0 & {\sin\quad é} & 0 & \quad & \quad & \quad & \quad & {x^{\prime} = {{\left( {\cos\quad é} \right)*x} + {\left( {\sin\quad é} \right)*z}}} \\0 & 1 & 0 & 0 & \quad & \quad & \quad & \quad & {y^{\prime} = y} \\{{- \sin}\quad é} & 0 & {\cos\quad é} & 0 & \quad & \quad & \quad & \quad & {z^{\prime} = {{{- \left( {\sin\quad é} \right)}*x} + {\left( {\cos\quad é} \right)*z}}} \\0 & 0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$

And the matrix for rotating around the z axis by angle é is:$\begin{matrix}{\cos\quad é} & {{- \sin}\quad é} & 0 & 0 & \quad & \quad & \quad & \quad & {x^{\prime} = {{\left( {\cos\quad é} \right)*x} - {\left( {\sin\quad é} \right)*y}}} \\{\sin\quad é} & {\cos\quad é} & 0 & 0 & \quad & \quad & \quad & \quad & {y^{\prime} = {{\left( {\sin\quad é} \right)*x} + {\left( {\cos\quad é} \right)*y}}} \\0 & 0 & 1 & 0 & \quad & \quad & \quad & \quad & {z^{\prime} = z} \\0 & 0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$Mirroring

Mirroring involves flipping all points about an arbitrary plane in 3Dspace, as illustrated in FIG. 15. From the diagram in FIG. 15 we caneasily see that to mirror a point we need to:

-   -   1) Calculate the distance of the point from the plane    -   2) Move the point in the opposite direction to the plane normal        a distance of 2*dist.

The distance of a point from a plane can be calculated easily with theplane equation:dist=p*normal+k

(Normally this equation assumes that the plane normal is a unit vector.In this particular case however this is not a requirement).

Based on the information above we can see that the final mirror equationis: $\begin{matrix}{p^{\prime} = {p - {2*{dist}*{normal}}}} \\{= {p - {2*\left( {{p*{normal}} + k} \right)*{normal}}}} \\{= {p - {2*\left( {{{p.x}*{{normal}.x}} + {{p.y}*{{normal}.y}} + {{p.z}*{{normal}.z}} + k} \right)*{normal}}}}\end{matrix}$

Expanding this out we get the equation for each element in p′:$\begin{matrix}{{p^{\prime} \cdot x} = {{p \cdot x} - {2*{p \cdot x}*{nx}*{nx}} + {2*{p \cdot y}*{ny}*{nx}} +}} \\{{2*{p \cdot z}*{nz}*{nx}} + {2*k*{nx}}} \\{{p^{\prime} \cdot y} = {{p \cdot y} - {2*{p \cdot x}*{nx}*{ny}} + {2*{p \cdot y}*{ny}*{ny}} +}} \\{{2*{p \cdot z}*{nz}*{ny}} + {2*k*{ny}}} \\{{p^{\prime} \cdot z} = {{p \cdot z} - {2*{p \cdot x}*{nx}*{nz}} + {2*{p \cdot y}*{ny}*{nz}} +}} \\{{2*{p \cdot z}*{nz}*{nz}} + {2*k*{nz}}}\end{matrix}$where <nx,ny,nz>=normal. Thus the final mirror matrix for any givenplane p*<nx,ny,nz>+k=0 is: $\begin{matrix}{1 - {2*{nx}*{nx}}} & {{- 2}*{nx}*{ny}} & {{- 2}*{nx}*{nz}} & {{- 2}*{nx}*k} \\{{- 2}*{ny}*{nx}} & {1 - {2*{ny}*{ny}}} & {{- 2}*{ny}*{nz}} & {{- 2}*{ny}*k} \\{{- 2}*{nz}*{nx}} & {{- 2}*{nz}*{ny}} & {1 - {2*{nz}*{nz}}} & {{- 2}*{nz}*k} \\0 & 0 & 0 & 1\end{matrix}$(Note  the  common  terms : m[0][1] = m[1][0], m[0][2] = m[2][0], and  m[1][2] = m[2][1]).Multiplying Matrices

So now that we know how to represent each of the different types ofmatrix operations, how do we combine them? By multiplying two matricestogether:P=B×A

If transforming a point by A gives one effect, and transforming it by Bgives another, then transforming it by P alone gives you the same resultas if you had transformed it by A then by B. You can keep multiplying amatrix by as many other matrices as you like, and each time the endproduct will contain the info for all of them in the correct order.

To multiply B by A, you treat each column in A as a separate columnvector, and transform it by matrix B to get the new column. Let's lookat doing it for the first column in A: ${\begin{matrix}{P11} \\{P21} \\{P31} \\{P41}\end{matrix}} = {{\begin{matrix}{B11} & {B12} & {B13} & {B14} \\{B21} & {B22} & {B23} & {B24} \\{B31} & {B32} & {B33} & {B34} \\{B41} & {B42} & {B43} & {B44}\end{matrix}}{\begin{matrix}{A11} \\{A21} \\{A31} \\{A41}\end{matrix}}}$

So that first column of P will be: $\quad{\begin{matrix}{\left( {{A11}*{B11}} \right) + \left( {{A21}*{B12}} \right) + \left( {{A31}*{B13}} \right) + \left( {{A41}*{B14}} \right)} \\{\left( {{A11}*{B21}} \right) + \left( {{A21}*{B22}} \right) + \left( {{A31}*{B23}} \right) + \left( {{A41}*{B24}} \right)} \\{\left( {{A11}*{B31}} \right) + \left( {{A21}*{B32}} \right) + \left( {{A31}*{B33}} \right) + \left( {{A41}*{B34}} \right)} \\{\left( {{A11}*{B41}} \right) + \left( {{A21}*{B42}} \right) + \left( {{A31}*{B43}} \right) + \left( {{A41}*{B44}} \right)}\end{matrix}}$

We need to split the A matrix up into it's 4 column vectors andtransform each column vector by matrix B to get 4 new column vectors:${\begin{matrix}{P11} \\{P21} \\{P31} \\{P41} \\{P12} \\{P22} \\{P32} \\{P42} \\{P13} \\{P23} \\{P33} \\{P43} \\{P14} \\{P24} \\{P34} \\{P44}\end{matrix}}\begin{matrix}\quad \\ = \\\quad \\\quad \\\quad \\ = \\\quad \\\quad \\\quad \\ = \\\quad \\\quad \\\quad \\ = \\\quad \\\quad\end{matrix}\quad{\begin{matrix}{B11} & {B12} & {B13} & {B14} \\{B21} & {B22} & {B23} & {B24} \\{B31} & {B32} & {B33} & {B34} \\{B41} & {B42} & {B43} & {B44} \\{B11} & {B12} & {B13} & {B14} \\{B21} & {B22} & {B23} & {B24} \\{B31} & {B32} & {B33} & {B34} \\{B41} & {B42} & {B43} & {B44} \\{B11} & {B12} & {B13} & {B14} \\{B21} & {B22} & {B23} & {B24} \\{B31} & {B32} & {B33} & {B34} \\{B41} & {B42} & {B43} & {B44} \\{B11} & {B12} & {B13} & {B14} \\{B21} & {B22} & {B23} & {B24} \\{B31} & {B32} & {B33} & {B34} \\{B41} & {B42} & {B43} & {B44}\end{matrix}}{\begin{matrix}{A11} \\{A21} \\{A31} \\{A41} \\{A12} \\{A22} \\{A32} \\{A42} \\{A13} \\{A23} \\{A33} \\{A43} \\{A14} \\{A24} \\{A34} \\{A44}\end{matrix}}$

The resulting matrix is then made by combing these 4 column vectors backinto a single matrix: ${B \times A} = {\begin{matrix}{P11} & {P12} & {P13} & {P14} \\{P21} & {P22} & {P23} & {P24} \\{P31} & {P32} & {P33} & {P34} \\{P41} & {P42} & {P43} & {P44}\end{matrix}}$

One thing to keep in mind is that the bottom row in any matrix willalways be [0 0 0 1], so we can use this to slightly speed up matrixmultiplication. Here's a general algorithm which will multiply 2matrices:

Let A and B be our matrices, P will be the result (B x A). var i,,j :integer; for i := 0 to 2 do  for j := 0 to 3 do   begin    P[i][j] := 0;   for k := 0 to 3 do     P[i][j] := P[i][j] + B[i][k] * A[k][j];   end;P[3][0] := 0; { Set the bottom row to 0 0 0 1 } P[3][1] := 0; P[3][2] :=0; P[3][3] := 1.The Inverse Matrix

Let's say we have two matrices A and B. If the following is true:${A \times B} = {{B \times A} = {\begin{matrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{matrix}}}$then we say that B is the _inverse_ of A (and visa-versa). If youtransform a point P1 by matrix A then you'll get a new point P2. If youthen transform P2 by matrix B, it'll return P1. The two matrixes havethe same effect but in “opposite” directions, e.g., if A moves allpoints 5 spaces to the “left” then B will move them 5 spaces to the“right”. Similarly if A rotates space around a particular axis one waythen B will rotate by the same amount around the same axis in theopposite direction.

There are several methods for calculating the inverse of a matrix, forlarge matrices the Gaussian method is preferred. Gaussian uses atechnique called of “row reduction”, first we create a large matrix likeso: $\begin{matrix}{{This}\quad{is}\quad{the}\quad{matrix}} & \quad \\{{we}\quad{are}\quad{trying}\quad{to}} & {{This}\quad{is}\quad{the}\quad{identity}} \\{{{find}\quad{the}\quad{inverse}\quad{of}}\quad} & {\quad{matrix}} \\{\begin{matrix}{A11} & {A12} & {A13} & {A14} \\{A21} & {A22} & {A23} & {A24} \\{A31} & {A32} & {A33} & {A34} \\{A41} & {A42} & {A43} & {A44}\end{matrix}} & {\quad\begin{matrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{matrix}}\end{matrix}$

Our goal is to perform various “elementary row operations” on thismatrix to put it in the following form: $\begin{matrix}{\quad{{Identity}\quad{matrix}}} & {{The}\quad{inverse}} \\{\begin{matrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{matrix}} & {\quad\begin{matrix}{B11} & {B12} & {B13} & {B14} \\{B21} & {B22} & {B23} & {B24} \\{B31} & {B32} & {B33} & {B34} \\{B41} & {B42} & {B43} & {B44}\end{matrix}}\end{matrix}$

In the above matrices A is our initial matrix and B will be the inverseof A.

There are three kinds of elementary row operations we can use:

-   -   RULE 1—Interchange two rows.    -   RULE 2 Multiply any row by a constant other than 0.    -   RULE 3—Add a constant multiple of any row to another row.

Here is the basic algorithm to figure out which steps to perform as wellas their order (this algorithm assumes that the bottom row of a matrixis always [0 0 0 1]):

var i,j,row : integer; multiple, divisor : real;  A : array[0..3][0..7]of real; Set values A[0..3][0..3] to our 4 x 4 matrix Set valuesA[0..3][4..7] to the 4 x 4 identity matrix { Loop through each row } forrow := 0 to 3 do  begin   { Make sure this row doesn't have a 0 inA[row][row] (use    RULE 1) }   if A[row][row] = 0 then    begin    find a row i where (i>row) and (i<3) and (A[i][row] <> 0)    interchange rows ‘i’ and ‘row’    end;   { Divide this row by aconstant so that A[row][row] = 1 (use RULE 2) }   divisor :=A[row][row];   for j := 0 to 7 do    A[row][j] := A[row][j] / divisor;  { Make all other elements in column ‘row’ a 0 (use RULE 3) }   for i:= 0 to 2 do    if i <> row then     begin      multiple := A[i][row];     for j := 0 to 7 do       A[i][j] := A[i][j] − multiple * A[row][j];    end;  end; Return the matrix in the values A[0..3][4..7]

There are cases where a matrix doesn't have an inverse. If this happensthen the first part of the algorithm above will fail, you won't be ableto find a row where A[i][row]!=0. However, if you stick to the “normal”3D operations discussed in this article (translation, rotation, scalingetc) then the matrices produced will always have an inverse.

One other very important thing to remember is that due to round-offerrors you will often get values in the matrix such as 0.000001 whichare supposed to be 0. When looking for a row where A[i][row]!=0 you mustkeep this in mind, and instead check that the absolute value is largerthan some small number (e.g., 0.0001).

Rotating Around an Arbitrary Axis

Let's say we need to rotate space by about an arbitrary axis defined bythe line o+kd (where o is the origin point <ox,oy,oz> and d is thedirectional vector <dx,dy,dz>). Here are the steps needed to performthis rotation:

1) First translate all of space so that the line of rotation starts atthe origin. Do this by subtracting the line's origin point from allpoints. Here is the matrix needed, along with it's inverse:$\begin{matrix}{F = {\begin{matrix}1 & 0 & 0 & {- {ox}} \\0 & 1 & 0 & {- {oy}} \\0 & 0 & 1 & {- {oz}} \\0 & 0 & 0 & 1\end{matrix}}} & {F^{\prime} = {\begin{matrix}1 & 0 & 0 & {ox} \\0 & 1 & 0 & {oy} \\0 & 0 & 1 & {o\quad z} \\0 & 0 & 0 & 1\end{matrix}}}\end{matrix}$

2) Next rotate space around the z axis until the line of rotation liesin the x/z plane: $\begin{matrix}{G = {\begin{matrix}{{dx}/v} & {{dy}/v} & 0 & 0 \\{{- {dy}}/v} & {{dx}/v} & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{matrix}}} & \quad & \quad & \quad & {G^{\prime} = {\begin{matrix}{{dx}/v} & {{- {dy}}/v} & 0 & 0 \\{{dy}/v} & {{dx}/v} & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{matrix}}}\end{matrix}$where v=SQRT(dx*dx+dy*dy). You will now have a new line of rotation withone end at the origin and the other end at the point <v,0,dz>.

3) Now rotate space around the y axis until the line of rotation liesalong the z axis: $\begin{matrix}{H = {\begin{matrix}{{dz}/w} & 0 & {{- v}/w} & 0 \\0 & 1 & 0 & 0 \\{v/w} & 0 & {{dz}/w} & 0 \\0 & 0 & 0 & 1\end{matrix}}} & \quad & \quad & \quad & {H^{\prime} = {\begin{matrix}{{dz}/w} & 0 & {v/w} & 0 \\0 & 1 & 0 & 0 \\{{- v}/w} & 0 & {{dz}/w} & 0 \\0 & 0 & 0 & 1\end{matrix}}}\end{matrix}$where w=SQRT(v*v+dz*dz)=SQRT(dx*dx+dy*dy+dz*dz).

4) At this point the axis of rotation is a line lying along the z axisbetween the points <0,0,0> and <0,0,w>, so you can rotate space aroundthe z axis by: $W = {\begin{matrix}{\cos\quad A} & {{- \sin}\quad A} & 0 & 0 \\{\sin\quad A} & {\cos\quad A} & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{matrix}}$

5) So to do the entire rotation you must first transform space by F,Gand H to align the axis of rotation along the z axis, then you do theactual rotation around the z axis using matrix W, then you transformspace by H′, G′ and F′ to put everything back into place. By multiplyingthese matrices together you can create a single matrix to do the entirerotation:P=F′×G′×H′×W×H×G×F

Any points transformed by this matrix will be rotated about the line byan angle of A.

If the term 1/v in the G and G′ matrices result in a singularity, it canbe fixed by replacing G and G′ with the identity matrix when v is 0.

Solving Simultaneous Equations

Imagine that we have the following equations: $\begin{matrix}{{{a*x} + {b*y} + {c*z} + d} = 0} \\{{{e*x} + {f*y} + {g*z} + h} = 0} \\{{{i*x} + {j*y} + {k*z} + l} = 0}\end{matrix}$

Now imagine that we need to find a point <x,y,z> which satisfies all 4equations. These equations contain 12 unknowns, so it's clear that itwould be very difficult to solve using “regular” algebra.

Matrices on the other hand provide us with a very useful tool forsolving such problems. We can create a matrix like this:${\begin{matrix}0 \\0 \\0 \\1\end{matrix}} = {{\begin{matrix}a & b & c & d \\e & f & g & h \\i & j & k & l \\0 & 0 & 0 & 1\end{matrix}}\quad{\begin{matrix}x \\y \\z \\1\end{matrix}}}$

When we expand this out (as I showed earlier on in this article) we getthe 4 equations we are trying to solve. We also know that we can use theinverse matrix to rearrange the equation: ${\begin{matrix}x \\y \\z \\1\end{matrix}} = {{\begin{matrix}a & b & c & d \\e & f & g & h \\i & j & k & l \\0 & 0 & 0 & 1\end{matrix}}^{- 1}\quad{\begin{matrix}0 \\0 \\0 \\1\end{matrix}}}$

These equations will expand out to the following: $\begin{matrix}{x = {{{A*0} + {B*0} + {C*0} + {D*1}} = D}} \\{y = {{{E*0} + {F*0} + {G*0} + {H*1}} = H}} \\{z = {{{I*0} + {J*0} + {K*0} + {L*1}} = L}}\end{matrix}$

(The letters have been capitalized since the inverse matrix willtypically contain different values than the original).

In summary, if we create a matrix from 3 parametric equations and thencalculate the inverse of that matrix then the 4th column will contain apoint which satisfies all 3 equations. If the inverse matrix operationfails then we know that the problem has no solution (or an infinitenumber of solutions).

Incidentally the fact that we'll only wind up using the 4th column canhelp in optimizing this procedure since we don't need to calculate theentire matrix.

2D Matrices

In this text I've been using 4×4 matrices for 3D coordinate geometry,but it's just as easy to use a smaller matrix for 2D (or higher), wesimply drop the z term.

2D matrices are of the form: ${\begin{matrix}x^{\prime} \\y^{\prime} \\1\end{matrix}} = {{\begin{matrix}{M11} & {M12} & {M13} \\{M21} & {M22} & {M23} \\0 & 0 & 1\end{matrix}}\quad{\begin{matrix}x \\y \\1\end{matrix}}}$ i.e., x^(′) = M11 * x + M12 * y + M13  y^(′) = M21 * x + M22 * y + M23

The 2identify matrix is: $\begin{matrix}1 & 0 & 0 & \quad & \quad & \quad & \quad & \quad & {x^{\prime} = x} \\0 & 1 & 0 & \quad & \quad & \quad & \quad & \quad & {y^{\prime} = y} \\0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$

To translate by <tx, ty>: $\begin{matrix}1 & 0 & {tx} & \quad & \quad & \quad & \quad & \quad & {x^{\prime} = {x + {tx}}} \\0 & 1 & {ty} & \quad & \quad & \quad & \quad & \quad & {y^{\prime} = {y + {ty}}} \\0 & 0 & 1 & \quad & \quad & \quad & \quad & \quad & \quad\end{matrix}\quad$

To scale by <sx, sy>: $\begin{matrix}{\quad\begin{matrix}{\begin{matrix}{sx} & 0 & 0 \\0 & {sy} & 0 \\0 & 0 & 1\end{matrix}} & \begin{matrix}{x^{\prime} = {x*{sx}}} \\{y^{\prime} = {y*{sy}}}\end{matrix}\end{matrix}\quad} & \quad\end{matrix}$

And to rotate about the origin by A radians: $\begin{matrix}{\begin{matrix}{\cos\quad(A)} & {{- \sin}\quad(A)} & 0 \\{\sin\quad(A)} & {\cos\quad(A)} & 0 \\0 & 0 & 1\end{matrix}} & {{x^{\prime} = {{\cos\quad(A)*x} - {\sin\quad(A)*y}}}{y^{\prime} = {{\sin\quad(A)*x} + {\cos\quad(A)*y}}}}\end{matrix}$

Copyright (c) 1997 Mark Feldman (pcgpe@geocities.com)—All RightsReserved

This article is part of The Win95 Game Programmer's Encyclooedia

Please retain this footer if you distribute this file.

Although the present invention and its advantages have been described indetail, it should be understood that various changes, substitutions andalterations can be made herein without departing from the spirit andscope of the invention. Moreover, the scope of the present applicationis not intended to be limited to the particular embodiments of theprocess, machine, manufacture, composition of matter, means, methods andsteps described in the specification. As one of ordinary skill in theart will readily appreciate from the disclosure of the presentinvention, processes, machines, manufacture, compositions of matter,means, methods, or steps, presently existing or later to be developedthat perform substantially the same function or achieve substantiallythe same result as the corresponding embodiments described herein may beutilized according to the present invention. Accordingly, the inventionis intended to include within its scope such processes, machines,manufacture, compositions of matter, means, methods, or steps.

1. In a digitized tomosynthesis method for projecting a 3D volumetricimage of an object onto a virtual projection plane, the image beingdefined by object voxels, in which a ray of energy from a source travelsthrough the object to impinge on an energy sensing image plate having anactive area defining an image plane and in which an image is acquired bythe energy sensor at successive relative rotational positions of theobject and source, wherein a Z axis is defined as the axis of relativerotation of the object, an X axis is defined as relatively horizontal atright angles to the Z axis, a Y axis is defined as relatively verticalat right angles to the X and Z axes, rotation about the X axis isdefined as pitch, rotation about the Y axis is defined as yaw, androtation about the Z axis is defined as roll, and wherein the imageplate is parallel with the XZ plane and the source is positioned abovethe XZ plane with its optical axis defining the center of a radiationcone and intersecting the image plate at approximately its centralpixel, the improvement according to which a 3D volumetric image of theobject is projected by steps comprising: computing the coordinate of anobject voxel referenced to the image plane; computing the image planeintercept of the ray of energy from the source through the object voxel;and projecting a reconstructed 3D volumetric image of the object ontothe virtual projection plane.
 2. The method of claim 1 in which thecoordinate of the object voxel referenced to the image plane is computedby: defining an object bounding box above the XZ plane that interceptsthe radiation cone so that a shadow of the object falls on the activearea of the image plate, the position and orientation of the boundingbox having a reference point coordinate defined by X, Y, Z, pitch, yaw,and roll; specifying the coordinate of the voxel referenced to theobject bounding box; rotating the voxel coordinate by the pitch, yaw,and roll of the bounding box; and translating the voxel coordinate tothe reference point coordinate.
 3. The method of claim 2 in which thebounding box has a size and height whereby it includes the volume ofinterest of the object.
 4. The method of claim 2 in which the boundingbox and image plate are each rectangular and in which the sides of thebounding box are parallel with the sides of the image plate.
 5. Themethod of claim 1 in which the virtual projection plane comprises analpha blend plane subdivided into an array of pixel cells, wherein theray of energy defines a trace thereof representing a virtual sight pathray passing from a view point through one pixel cell of the alpha blendplane to define a view level image.
 6. The method of claim 5 in whichthe array of pixel cells is a rectangular array.
 7. The method of claim5 in which there is a ray trace for each pixel cell of the alpha blendplane for each view level image.
 8. The method of claim 5 comprisingforming a projected planar view of a selected view level image plane atthe alpha blend plane by transferring the image pixel cell value at theintersect of the ray trace and the selected view level plane to theimage pixel cell at the intersect of the ray trace and the alpha blendplane.
 9. The method of claim 1 in which the steps thereof are performedmathematically and programmed whereby the reconstructed 3D volumetricimage is projected by execution of the program on a computer.
 10. Themethod of claim 1 in which a plurality of view level images are formedsuccessively at one image plane.
 11. The method of claim 10 wherein theview point and alpha blend plane are translated a level closer to saidone image plane for each successive view level.
 12. The method of claim1 in which the 3D volumetric image is stored in a computer memory as anarray of pixels represented by a value.
 13. The method of claim 12 inwhich the brightness of each pixel is proportional to the value.
 14. Themethod of claim 1 in which the ray of energy is x-ray radiation.
 15. Themethod of claim 1 in which the energy sensor is a flat panel digitaldetector.
 16. In a digitized tomosynthesis method for projecting a 3Dvolumetric image of an object onto a virtual projection plane, the imagebeing defined by object voxels, in which a ray of energy from a sourcetravels through the object to impinge on a rectangular energy sensingimage plate having an active area defining an image plane and in whichan image is acquired by the energy sensor at successive relativerotational positions of the object and source, wherein a Z axis isdefined as the axis of relative rotation of the object, an X axis isdefined as relatively horizontal at right angles to the Z axis, a Y axisis defined as relatively vertical at right angles to the X and Z axes,rotation about the X axis is defined as pitch, rotation about the Y axisis defined as yaw, and rotation about the Z axis is defined as roll, andwherein the image plate is parallel with the XZ plane and the source ispositioned above the XZ plane with its optical axis defining the centerof a radiation cone and intersecting the image plate at approximatelyits central pixel, the improvement according to which the virtualprojection plane comprises an alpha blend plane subdivided into arectangular array of pixel cells, wherein the ray of energy defines atrace thereof representing a virtual sight path ray passing from a viewpoint through one pixel cell of the alpha blend plane to define a viewlevel image, there being a ray trace for each pixel cell for each viewlevel image, and in which a 3D volumetric image of the object isprojected by steps comprising: computing the coordinate of an objectvoxel referenced to the image plane by: (a)  defining a rectangularobject bounding box above the XZ plane, the bounding box having sidesthat are parallel with the sides of the image plate and having a sizeand height that includes the volume of interest of the object, wherebythe bounding box intercepts the radiation cone so that a shadow of theobject falls on the active area of the image plate, the position andorientation of the bounding box having a reference point coordinatedefined by X, Y, Z, pitch, yaw, and roll, (b) specifying the coordinateof the voxel referenced to the object bounding box, (c) rotating thevoxel coordinate by the pitch, yaw, and roll of the bounding box, and(d) translating the voxel coordinate to the reference point coordinate;computing the image plane intercept of the ray of energy from the sourcethrough the object voxel; and forming a plurality of successivelyselected view level image planes at the alpha blend plane bytransferring the image pixel cell values at the intersects of the raytrace and the selected view level planes to the image pixel cells at theintersects of the ray trace and the alpha blend plane, the view pointand alpha blend plane being translated a level closer to said one imageplane for each successive view level, whereby to project a reconstructed3D volumetric image of the object onto the virtual projection plane. 17.The method of claim 16 in which the steps thereof are performedmathematically and programmed whereby the reconstructed 3D volumetricimage is projected by execution of the program on a computer.
 18. Themethod of claim 16 in which the 3D volumetric image is stored in acomputer memory as an array of pixels represented by a value in whichthe brightness of each pixel is proportional to the value.
 19. Themethod of claim 16 in which the ray of energy is x-ray radiation. 20.The method of claim 16 in which the energy sensor is a flat paneldigital detector.